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Abstract 

We present a numerical procedure of solving the subdiffusion equation 
with Caputo fractional time derivative. On the basis of few examples we 
show that the subdiffusion is a 'long time memory' process and the short 
memory principle should not be used in this case. 

PACS numbers: 02.50.Ey, 05.10.-a, 02.60.Cb 

1 Introduction 

The subdiffusion equation is of fractional order with respect to the time vari- 
able. Unfortunately, the exact solutions are known only for relatively simple 
systems, similarly to the normal diffusion case. In more complicated situations 
such as a system with subdiffusion coefficient depending on concentration, in- 
homogenous fractional subdiffusive system or subdiffusion-reaction system, one 
needs a numerical procedure to solve of the equation. Usually the subdiffusion 
equation contains the Riemann-Liouville fractional time derivative of the order 
1 — a (a denotes here the subdiffusion parameter), which is not convenient for 
physical interpretation of initial conditions. To get the subdiffusion equation 
with initial conditions which have simple interpretation, one can use the subd- 
iffusion equation with the Caputo fractional time derivative of the order a. As 
far as we know, there is a numerical method to solve the subdiffusion equation 
with the Riemann-Liouville fractional time derivative The equation with 
Caputo derivative has been numerically studied only within the time fractional 
discrete random walk (2). 

Subdiffusion is a process with the time memory. There arises a practical 
problem with the memory length, which extends to — oo. To omit the difficulty 
Podlubny 01 |S] postulated to apply the short memory principle which assumes 
that the relatively small memory length is sufficient to obtain satisfactorily 
accuracy of numerical solutions for sufficiently long times. However, as shown 
here, the method produces significant differences between the numerical and 



exact solutions. In this paper we present a numerical procedure of solving the 
subdifFusion equation with Caputo derivative, which is based on the fractional 
difference approach, and we briefly study the efficiency of the short memory 
principle for this case. 



2 Subdiffusion equation 

The transport process is described by the subdiffusion equation ^ 

dC{x,t) _ "^8^-°' d'^C{x,t) 

dt ^ " dx^ ' ^ ' 

where < a < 1, C{x,t) denotes the concentration of transported substance, 
the Riemann-Liouville fractional time derivative is defined for a > as 

^^^"/w ^ 1 r\,. fit') 

T{n-a)dt''Jo 

where the integer number n fulfills the relation n — 1 < a < n. 

The presence of time derivatives on both sides of Eq. Q is not convenient 
for numerical calculations. To simplify the numerical procedure we rewrite the 
Eq. (P) in the form 

_ n d^C{x,t) 

where the fractional time derivative on the left-hand side of Eq. ^ is now 
Caputo derivative defined by the relation 

dt- T{a-n)Jo 

/("^ denotes the derivative of natural order n. 

We note that the Laplace transform of the Caputo fractional derivative in- 
volves values of the function / and its derivatives of natural order at i = while 
Laplace transform of the Riemann-Liouville fractional derivative includes the 
fractional derivatives of / at t = 0. Thus, the physical interpretation of initial 
conditions in the former case is clear in contrast to the latter one So, it 
is more convenient to set the initial conditions for the equation with Caputo 
derivative than for the equation with Riemann-Liouville one. 



3 Numerical procedure 
3.1 Fractional derivatives 

To numerically solve the normal diffusion equation one usually substitutes the 

time derivative by the backward difference — '^''*^^At "'^^^ • ^^'^ 

sented procedure we proceed in a similar way. For that purpose we use the 
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Griinwald-Letnikow fractional derivative, which is defined as a hmit of a fractional- 
order backward difference ^ 

^^-i!";o<^')-°E<-i)'(;)/('-'"^'). (5) 

r=0 ^ ^ 

where a > 0, [x] means the integer part of x and 



a 



T{a + 1) _ a{a - l){a - 2) • . . . • [a - (r - 1)] 



r J rW{a - r + l) l-2-3-...-r 

When the function f{t) of positive argument has continuous derivatives of 
the integer order 0, 1, . . . ,n, the Riemann-Liouville definition (j^J is equivalent 
to the Griinwald-Letnikow one 0]. So, we can take 

The relation between Riemann-Liouville and Caputo derivatives is more com- 
plicated and reads as 

— — = ^^ + Z^*fc-"+iw/ (0)' (") 

where 

i>0 



From Eqs. ©-((HI) we can express the Caputo fractional derivative in terms of 
the fractional-order backward difference ^ 



at" At^o' ' V ^ / <"r(i-Q;) 



r=0 

3.2 Algorithm 

The standard way to approximate of the fractional derivative, which is useful 
for numerical calculations, is to omit the limit in Eq. lO and to change the 
infinite series occurring in © to the finite one 

^ - (Ar- h-^r ( : ) /(^ - r^i) - j^^^m, m 

r=0 ^ ^ ^ ' 

where arbitrary chosen parameter L is called the memory length. Substituting 
Eq. (|10|l to Eq. and using the following approximation of the second order 
derivative 

d^fjx) fix + Ax) - 2 fix) + fix - Ax) 
dx^ ~ (Aa;)2 ' ^ 
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after simple calculation we obtain 



r—l 

1 (At)" 

C{x,0) + Daj-^[C{x + Ax,t~At) (12) 



i"r(l-a)^ "(Aa;)2 

-2C(a;, f - At) + C(a: - Ax, i - At)]. 

There arise a problem with the choice of initial conditions. We choose the 
initial conditions only for t = assuming that the concentration given for earlier 
moments t < does not influence the process for t > 0. This assumption is in 
agreement with the procedure of solving the equations with Caputo fractional 
derivative where the initial conditions are determined only at t = for the 
derivatives of natural order f^^\t)\t=o, n = 0, 1, . . . , [a] (here /(°^ = /). In our 
considerations we have < a < 1, hence it is enough to set C{x, 0) as the initial 
condition. 

Starting with the initial condition C{x,0) we will find the time iterations 
C{x,tsAt) for ts = 1,2, ...,ts,max- When the number of time steps ts is less 
then the memory length L then we put L = ts in the series occurring in Eq. 
H12|) . otherwise the memory length is equal to L. 



4 Numerical results 



To test the numerical procedure we are going to compare the numerical solutions 
of the subdiffusion equation with the exact analytical ones. For that purpose 
we choose the homogenous system with the initial concentration 



C{x,0) = 



Co X <0 
a; > 



(13) 



The solution of the subdiffusion equation © with the initial condition H1,S|I is 
following [B] 



r< Ca ttW 



C{x,t) = < 



( [-^f- 
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X < 



X > 



(14) 



where H denotes the Fox function, which can be expressed by the series [7| 

2/a 1 1 \ o oo 1 / „ \ fe 



Dll'-t 



1 1 

2/a 



1 



a^^k\T{l-ka/2) 



(15) 



The results of numerical calculations and the analytical solutions are shown 
in the plots. In Figures 1-3 we present the numerical solutions of the subdiffusion 
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Figure 1: The concentration profiles C{x,t) calculated for a = 0.5, Da — 0.25, 
ts,max = 100 and with different memory length L = 100(0), 80(A), 50(o), 
30(V), and 10; continuous line represents the exact analytical solution. 

equation for different values of t, a, and Da- In each case we present the plot 
of analytical solution (continuous line) and numerical solutions calculated for 
different memory length L (symbols without line). The time and the memory 
length are given as the number of all time steps ts^max, which corresponds to 
the 'real time' t by the relation t = At ■ ts ^ax- In all cases we take Cq — 1, 
At — 0.1 and Ax ~ 0.5 (all quantities are given in arbitrary units); to calculate 
the analytical solutions H14|l we took 100 first terms in the series occurring in 
H15|l . We can see that the memory length determines the accuracy of numerical 
solutions. 

5 Final remarks 

We have presented the procedure to numerically solve the subdiffusion equation 
with Caputo fractional time derivative. The choice of the equation in such a 
form is not accidental since the interpretation of the initial condition in this 
case is simpler than in the equation with Riemann-Liouvillc derivative. In all 
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Figure 2: The concentration profiles obtained for a — 0.2, — 0.2, ts^max = 50 
and with memory length L = 50(^), 40(A), 30(o), 20(V), 10, and 5(x). 

considered cases the numerical solutions coincide with the analytical ones. In 
the studies ^E] the 'short memory principle' was postulated. According to this 
principle, the fractional derivative is approximated by the fractional derivative 
with moving lower limit t — L, where L is the 'memory length'. The examples 
presented in (4j suggest that the L = 50 time steps gives a good approximation 
for times of the order of ts ^ 100 time steps. However, the results presented 
here show that this memory length is not sufficient for the subdiffusion case. 
Our analysis demonstrate that the memory length should be longer than about 
80 per cent of the value of time variable. 
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Figure 3: The concentration profiles calculated for a = 0.3, Da = 0.1, time 
ts max = 50 and with memory length L = 50(0), 40(A), 30(o), 20(V), 10, and 
5(x). 
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